Mushroom Data Analysis

Índice

1. Introducción
   
2. Cargar librerías
   
3. Cargar datos
   
4. Visualización y preprocesamiento
   4.1 Imputación de valores nulos
   4.2 Codificación de variables categóricas
   4.3 Escalado de variables numéricas
5. Análisis de datos
   5.1 Análisis aprendizaje supervisado
   5.2 Análisis aprendizaje no supervisado
   5.2.1 K-means
   5.2.2 Clustering jerárquico
   

Introducción

Este trabajo se centra en el análisis de un conjunto de datos que contiene información sobre 61069 champiñones diferentes. Los datos han sido obtenidos de Kaggle y cada una de las instancias incluye 21 variables que describen diferentes aspectos de los champiñones, como su forma, tamaño y hábitat. En estas variables se encuentra inlcuida la clase de cada champiñón, indicando si es venenoso o comestible.

Antes de comenzar el análisis, es necesario realizar una fase de preprocesamiento de los datos. Esta fase tiene como objetivo preparar los datos para su análisis y obtener mayor conocimiento sobre estos. Durante esta fase, realizaremos tareas como visualizar los datos, detectar y tratar valores faltantes, o normalizar los datos. Una vez preparados los datos, procederemos a realizar el análisis. Para ello, utilizaremos tanto técnicas de aprendizaje supervisado como no supervisado.

El aprendizaje supervisado implica entrenar un modelo con datos etiquetados, es decir, que ya conocemos la clase de cada instancia. Una vez entrenado el modelo, podemos utilizarlo para predecir la clase de nuevas instancias. En este caso, utilizaremos diferentes algoritmos de clasificación para evaluar diferentes modelos.

El aprendizaje no supervisado, por otro lado, no requiere de datos etiquetados. En este caso, el objetivo es agrupar las instancias en diferentes grupos de forma que las instancias de un mismo grupo sean similares entre sí y diferentes a las de los demás grupos. Utilizaremos dos algoritmos de clustering para evaluar diferentes modelos. Además, cabe destacar, que en nuestro caso contamos con un conocimiento a priori de los datos, ya que conocemos la clase de cada instancia. Por tanto, podemos utilizar este conocimiento para evaluar los modelos de clustering.

A continuación se detalla los posibles valores que pueden tomar las variables del dataset:

1. class: edible=e, poisonous=p
2. cap-diameter: float number in cm
3. cap-shape: bell=b, conical=c, convex=x, flat=f, knobbed=k, sunken=s
4. cap-surface: fibrous=f, grooves=g, scaly=y, smooth=s
5. cap-color: brown=n, buff=b, cinnamon=c, gray=g, green=r, pink=p, purple=u, red=e, white=w, yellow=y
6. does.bruise.or.bleed: bruises-or-bleeding=t,no=f
7. gill-attachment: attached=a, descending=d, free=f, notched=n
8. gill-spacing: close=c, crowded=w, distant=d
9. gill-color: black=k, brown=n, buff=b, chocolate=h, gray=g, green=r, orange=o, pink=p, purple=u, red=e, white=w, yellow=y
10. stem-height: float number in cm
11. stem-width: float number in mm
12. stem-root: bulbous=b, swollen=s, club=c, cup=u, equal=e, rhizomorphs=z, rooted=r
13. stem-surface: see cap-surface + none=f
14. stem-color: see cap-color + none=f
15. veil-type: partial=p, universal=u
16. veil-color: see cap-color + none=f
17. has-ring: ring=t, none=f
18. ring-type: cobwebby=c, evanescent=e, flaring=r, grooved=g, large=l, pendant=p, sheathing=s, zone=z, scaly=y, movable=m, none=f, unknown=?
19. spore-print-color: see cap color
20. habitat: grasses=g, leaves=l, meadows=m, paths=p, heaths=h, urban=u, waste=w, woods=d
21. season: spring=s, summer=u, autumn=a, winter=w

Autores

Este trabajo ha sido realizado por:

• Ana Díaz Muñoz

• María Isabel Ramos Blanco

• Deyan Rosenov Stanchev

• Javier Vilariño Mayo

Todos los integrantes del grupo han realizado de forma conjunta la visualización y el preprocesamiento de los datos. Ana y Deyan se ha encargado de realizar el aprendizaje supervisado, mientras que Javier y María Isabel se han encargado de realizar el aprendizaje no supervisado. Por otro lado, todos los componentes han investigado acerca del clustering utilizado en BigML.

Cargar librerías

En los siguientes fragmentos se incluye el código necesario para instalar y cargar las librerías necesarias para el análisis de los datos realizado.

#install.packages("caret")
#install.packages("tidyverse")
#install.packages("plotly")
#install.packages("dplyr")
#install.packages("factoextra")
#install.packages("dendextend")

library(caret)
library(tidyverse)
library(plotly)
library(dplyr)
library(cluster)
library(factoextra)
library(dendextend)

Cargar datos

En primer lugar, realizamos la lectura del fichero csv, separando las columnas por “;” y mostramos las primeras 6 filas del fichero. Este código nos permitirá acceder a los datos de champiñones y trabajar con ellos en nuestro código R.

mushroom <- read.csv("./data/data.csv", sep = ";")
head(mushroom)

Echamos un vistazo a las características de los atributos del dataset. En el caso de las variables numéricas, se puede observar valores como el mínimo, máximo, media, desviación estándar, etc. Por otro lado, en el caso de las variables categóricas no obtenemos información relevante.

summary(mushroom)

    class            cap.diameter     cap.shape         cap.surface         cap.color         does.bruise.or.bleed
 Length:61069       Min.   : 0.380   Length:61069       Length:61069       Length:61069       Length:61069        
 Class :character   1st Qu.: 3.480   Class :character   Class :character   Class :character   Class :character    
 Mode  :character   Median : 5.860   Mode  :character   Mode  :character   Mode  :character   Mode  :character    
                    Mean   : 6.734                                                                                
                    3rd Qu.: 8.540                                                                                
                    Max.   :62.340                                                                                
 gill.attachment    gill.spacing        gill.color         stem.height       stem.width      stem.root        
 Length:61069       Length:61069       Length:61069       Min.   : 0.000   Min.   :  0.00   Length:61069      
 Class :character   Class :character   Class :character   1st Qu.: 4.640   1st Qu.:  5.21   Class :character  
 Mode  :character   Mode  :character   Mode  :character   Median : 5.950   Median : 10.19   Mode  :character  
                                                          Mean   : 6.582   Mean   : 12.15                     
                                                          3rd Qu.: 7.740   3rd Qu.: 16.57                     
                                                          Max.   :33.920   Max.   :103.91                     
 stem.surface        stem.color         veil.type          veil.color          has.ring          ring.type        
 Length:61069       Length:61069       Length:61069       Length:61069       Length:61069       Length:61069      
 Class :character   Class :character   Class :character   Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character   Mode  :character   Mode  :character   Mode  :character  
                                                                                                                  
                                                                                                                  
                                                                                                                  
 spore.print.color    habitat             season         
 Length:61069       Length:61069       Length:61069      
 Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character  
                                                         
                                                         
                                                         

Visualización y preprocesamiento

En primer lugar, vamos a comprobar si existen valores nulos en el dataset. Para ello, utilizaremos la función colSums(is.na(mushroom)), que nos devolverá la suma de valores nulos de cada variable.

colSums(is.na(mushroom))

               class         cap.diameter            cap.shape          cap.surface            cap.color 
                   0                    0                    0                    0                    0 
does.bruise.or.bleed      gill.attachment         gill.spacing           gill.color          stem.height 
                   0                    0                    0                    0                    0 
          stem.width            stem.root         stem.surface           stem.color            veil.type 
                   0                    0                    0                    0                    0 
          veil.color             has.ring            ring.type    spore.print.color              habitat 
                   0                    0                    0                    0                    0 
              season 
                   0 

Según el resultado obtenido anteriormente, no existen valores nulos en el dataset. Sin embargo, si observamos el dataset, podemos observar que existen valores vacíos. Para poder trabajar con ellos, y que no nos de problemas a la hora de realizar el preprocesamiento, sustituiremos dichos valores vacíos por NA.

mushroom[mushroom == ""] <- NA

Comprobamos que se han sustituido correctamente los valores vacíos por NA, pudiendo comprobar la cantidad de valores nulos que hay en cada variable. Esta información nos ayudará a decidir si será conveniente eliminar dichas variables o no.

colSums(is.na(mushroom))

               class         cap.diameter            cap.shape          cap.surface            cap.color 
                   0                    0                    0                14120                    0 
does.bruise.or.bleed      gill.attachment         gill.spacing           gill.color          stem.height 
                   0                 9884                25063                    0                    0 
          stem.width            stem.root         stem.surface           stem.color            veil.type 
                   0                51538                38124                    0                57892 
          veil.color             has.ring            ring.type    spore.print.color              habitat 
               53656                    0                 2471                54715                    0 
              season 
                   0 

Podemos observar que existen 5 variables donde más del 50% de los valores son nulos. Estas son: stem.surface, veil.color, spore.print.color, stem.root, veil.type.

En este caso, decidimos eliminar aquellas variables que cuentan con más del 50% de sus valores faltantes. La razón de esta decisión es que, si decidimos imputar los valores faltantes, estaríamos inventando demasiados datos. Imputar valores faltantes significa reemplazar los valores faltantes con algún valor que estimemos adecuado. Sin embargo, si una variable tiene más del 50% de sus valores faltantes, significa que estaríamos reemplazando más de la mitad de los valores de esa variable. Esto nos llevaría a tener un conjunto de datos con demasiados valores inventados, lo que podría afectar la precisión de nuestro análisis. Para ello, en primer lugar, obtendremos el nombre de las columnas que queremos eliminar.

nacols <- colnames(mushroom)[colSums(is.na(mushroom)) > nrow(mushroom) / 2]
print(nacols)

[1] "stem.root"         "stem.surface"      "veil.type"         "veil.color"        "spore.print.color"

A continuación, eliminaremos dichas columnas del dataset.

mushroom <- mushroom[, !names(mushroom) %in% nacols]

Comprobamos que se han eliminado correctamente las columnas.

print(colnames(mushroom))

 [1] "class"                "cap.diameter"         "cap.shape"            "cap.surface"          "cap.color"           
 [6] "does.bruise.or.bleed" "gill.attachment"      "gill.spacing"         "gill.color"           "stem.height"         
[11] "stem.width"           "stem.color"           "has.ring"             "ring.type"            "habitat"             
[16] "season"              

Para poder analizar de forma específica cada variable, separamos las variables numéricas de las categóricas.

colsnames <- colnames(mushroom)
numerical_features <- c("cap.diameter", "stem.height", "stem.width")
categorical_features <- colsnames[!colsnames %in% numerical_features]
print(categorical_features)

 [1] "class"                "cap.shape"            "cap.surface"          "cap.color"            "does.bruise.or.bleed"
 [6] "gill.attachment"      "gill.spacing"         "gill.color"           "stem.color"           "has.ring"            
[11] "ring.type"            "habitat"              "season"              

print(numerical_features)

[1] "cap.diameter" "stem.height"  "stem.width"  

Comenzamos analizando las variables categóricas. Para ello visualizamos la distribución de las variables categóricas a través de histogramas. Observaremos los posibles valores de cada variable categórica junto con su frecuencia de aparición en el dataset.

for (i in categorical_features) {
  print(ggplot(mushroom, aes_string(x = i)) +
    geom_bar(fill = "#7fd6d9") +
    geom_text(stat = "count", aes(label = scales::percent(..count.. / nrow(mushroom)), vjust = -0.25)) +
    labs(x = i, y = "Percentage") +
    theme(axis.text.x = element_text(angle = 90, hjust = 1)))
}

Tras visualizar los histogramas de cada variable, obtenemos las siguientes conclusiones:

• La variable dependiente “class” está balanceada, es decir, la frecuencia de aparición de “e” (edible) y “p” (poisonous) es similar, en una proporción de 44,5% y 55,5% respectivamente. Por lo tanto, no es necesario realizar un balanceo de la variable dependiente, permitiendo aplicar la medida de evaluación “accuracy” para evaluar el modelo.
• Siguen existiendo valores NA en algunas variables, pero estos serán imputados posteriormente con la moda de cada variable.
• La variable ring-type tiene 8 posibles valores, pero en el 80% de los casos, el valor es “f”, por lo que quizás se podría eliminar esta variable en un futuro, ya que parece que aporta información irrelevante. Más adelante la analizaremos en profundidad.

A continuación, seguiremos con el análisis de las variables numéricas, donde visualizaremos su distribución respecto a la clase a través de una gráfica generada con “featurePlot”. Esta requiere que la variable objetivo sea de tipo factor, por lo que hacemos la conversión.

mushroom$class <- as.factor(mushroom$class)
featurePlot(x = mushroom[, numerical_features], y = mushroom$class, plot = "strip")

Con la gráfica anterior, se puede observar que para valores altos de cap.diameter, stem.height y stem.width, la probabilidad de que la clase sea “e” es mayor que la de “p”. Por lo tanto, podemos deducir que si una seta es de tamaño grande, su probabilidad de ser comestible es bastante mayor.

Anteriormente hemos comentado que la variable “ring.type” tiene 8 posibles valores, pero en casi un 80% de los casos, el valor es “f”. Para evaluar su posible eliminación o la de alguna otra variable, se puede utilizar la función “nearZeroVar”. Esta función devuelve un vector con los índices de las variables que tienen una varianza cercana a 0.

near_zero_col <- nearZeroVar(mushroom, saveMetrics = FALSE)
colnames(mushroom)[c(near_zero_col)]

[1] "ring.type"

Como suponíamos, la variable “ring.type” es la que tiene una varianza cercana a 0, y por tanto podría ser eliminada. A continuación, visualizaremos la correlación existente con la variable dependiente “class” para tomar una decisión.

print(ggplot(mushroom, aes_string(x = "ring.type")) +
  geom_bar(aes(fill = class)))

Tras visualizar la gráfica anterior, se puede observar que la distribución de la variable dependiente “class” es similar en casi todos los valores de “ring.type”. Sin embargo, en el caso de “ring.type” = “z” y “ring.type” = “m”, la distribución de la variable dependiente “class” es diferente. Por lo tanto, se decide finalemnte mantener la variable “ring.type”.

Imputación de valores nulos

Como se ha comentado anteriormente, los valores nulos de las variables categóricas se imputarán a través de la moda de cada variable. Este proceso lo hacemos de forma manual, ya que la función preProcess() de la librería caret no permite imputar valores nulos de variables categóricas.

for (i in categorical_features) {
  mushroom[, i][is.na(mushroom[, i])] <- names(which.max(table(mushroom[, i])))
}

Comprobamos que ya no existen valores nulos en las variables categóricas.

colSums(is.na(mushroom))

               class         cap.diameter            cap.shape          cap.surface            cap.color 
                   0                    0                    0                    0                    0 
does.bruise.or.bleed      gill.attachment         gill.spacing           gill.color          stem.height 
                   0                    0                    0                    0                    0 
          stem.width           stem.color             has.ring            ring.type              habitat 
                   0                    0                    0                    0                    0 
              season 
                   0 

Escalado de variables numéricas

Con el objetivo de que todas las variables tengan la misma escala y evitar que una variable tenga más peso que otra, se escalarán las variables numéricas. Para realizar este escalado, se utilizará la función preProcess() de la librería caret. Esta función devuelve un objeto de tipo “preProcess” que contiene la información necesaria para escalar las variables numéricas. A continuación, se realizará el escalado y se sustituirán las variables numéricas originales por las escaladas.

range_numeric <- preProcess(mushroom[, numerical_features], method = c("range"))
mushroom[, numerical_features] <- predict(range_numeric, newdata = mushroom[, numerical_features])
str(mushroom)

'data.frame':   61069 obs. of  16 variables:
 $ class               : Factor w/ 2 levels "e","p": 2 2 2 2 2 2 2 2 2 2 ...
 $ cap.diameter        : num  0.24 0.262 0.221 0.223 0.23 ...
 $ cap.shape           : chr  "x" "x" "x" "f" ...
 $ cap.surface         : chr  "g" "g" "g" "h" ...
 $ cap.color           : chr  "o" "o" "o" "e" ...
 $ does.bruise.or.bleed: chr  "f" "f" "f" "f" ...
 $ gill.attachment     : chr  "e" "e" "e" "e" ...
 $ gill.spacing        : chr  "c" "c" "c" "c" ...
 $ gill.color          : chr  "w" "w" "w" "w" ...
 $ stem.height         : num  0.5 0.53 0.525 0.465 0.487 ...
 $ stem.width          : num  0.164 0.175 0.171 0.154 0.166 ...
 $ stem.color          : chr  "w" "w" "w" "w" ...
 $ has.ring            : chr  "t" "t" "t" "t" ...
 $ ring.type           : chr  "g" "g" "g" "p" ...
 $ habitat             : chr  "d" "d" "d" "d" ...
 $ season              : chr  "w" "u" "w" "w" ...

Análisis de datos

Análisis aprendizaje supervisado

Para trabajar con el dataset mediante el aprendizaje supervisado; es decir, utilizando datos que son etiquetados mediante la intervención de un ser humano, utilizaremos diferentes tipos de clasificadores, algunos de ellos ya trabajados durante las clases prácticas de la asignatura y otros que eran desconocidos para nosotros y sobre los cuales hemos tenido que investigar anteriormente sobre su funcionamiento en R. Los algoritmos de clasificación que hemos seleccionado y los cuales vamos a aplicar son:

1. Regresión logística.

2. KNN o el Vecino más cercano.

3. Árboles de decisión.

4. Random Forest.

5. MSV o Máquina de Soporte Vectorial.

Una vez aplicados cada uno de los clasificadores, los compararemos entre sí y seleccionaremos el o los algoritmos que mayor precisión proporcionen sin llegar al sobreajuste, tratando de buscar que el resultado final sea generalizado para los datos.

En primer lugar, ante de comenzar a aplicar los clasificadores, dividiremos el dataset en dos conjuntos: el primer conjunto será el de entrenamiento o training y el segundo será el conjunto de prueba o test. El conjunto de entrenamiento lo utilizaremos para entrenar los distintos modelos y el conjunto de prueba lo utilizaremos para evaluar cada uno de ellos una vez obtenidos.

library(caTools)
set.seed(18)

split <- sample.split(mushroom$class, SplitRatio = 0.8)
training_set <- subset(mushroom, split == TRUE)
test_set <- subset(mushroom, split == FALSE)

table(training_set$class)

    e     p 
21745 27110 

table(test_set$class)

   e    p 
5436 6778 

Regresión Logística

La regresión logística permite predecir el resultado de una variable categórica en función de las variables independientes o predictoras. A continuación se muestra un resumen del conjunto de datos de entrenamiento con el cual vamos a trabajar una vez aplicada la función glm().

rl_classiffier <- glm(class ~ ., family = binomial, data = training_set)
summary(rl_classiffier)

Call:
glm(formula = class ~ ., family = binomial, data = training_set)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.98113  -0.77757   0.00045   0.75561   2.99676  

Coefficients: (2 not defined because of singularities)
                        Estimate Std. Error z value Pr(>|z|)    
(Intercept)           -16.146622 310.494270  -0.052 0.958526    
cap.diameter           -3.097299   0.273141 -11.340  < 2e-16 ***
cap.shapec             -1.587320   0.089853 -17.666  < 2e-16 ***
cap.shapef             -1.548574   0.056840 -27.244  < 2e-16 ***
cap.shapeo             -0.505473   0.102060  -4.953 7.32e-07 ***
cap.shapep             -1.355581   0.077525 -17.486  < 2e-16 ***
cap.shapes             -1.914206   0.067428 -28.389  < 2e-16 ***
cap.shapex             -1.587070   0.052730 -30.098  < 2e-16 ***
cap.surfacee            0.777609   0.081079   9.591  < 2e-16 ***
cap.surfaceg           -0.469069   0.067815  -6.917 4.62e-12 ***
cap.surfaceh           -0.869361   0.063758 -13.635  < 2e-16 ***
cap.surfacei            1.878605   0.116592  16.113  < 2e-16 ***
cap.surfacek            3.205155   0.106739  30.028  < 2e-16 ***
cap.surfacel           -1.338110   0.100947 -13.256  < 2e-16 ***
cap.surfaces           -0.913831   0.059362 -15.394  < 2e-16 ***
cap.surfacet           -0.037125   0.050461  -0.736 0.461905    
cap.surfacew           -0.483820   0.084640  -5.716 1.09e-08 ***
cap.surfacey           -0.579375   0.064113  -9.037  < 2e-16 ***
cap.colore              1.986990   0.109951  18.072  < 2e-16 ***
cap.colorg              0.582318   0.108490   5.367 7.98e-08 ***
cap.colork              1.553011   0.131690  11.793  < 2e-16 ***
cap.colorl             -0.129820   0.140094  -0.927 0.354103    
cap.colorn              0.374666   0.099619   3.761 0.000169 ***
cap.coloro              1.627975   0.112745  14.439  < 2e-16 ***
cap.colorp              1.051060   0.127472   8.245  < 2e-16 ***
cap.colorr              2.925006   0.136340  21.454  < 2e-16 ***
cap.coloru              1.538956   0.121348  12.682  < 2e-16 ***
cap.colorw              0.861936   0.104891   8.217  < 2e-16 ***
cap.colory              0.615058   0.102936   5.975 2.30e-09 ***
does.bruise.or.bleedt  -0.147991   0.038250  -3.869 0.000109 ***
gill.attachmentd        0.644857   0.046979  13.726  < 2e-16 ***
gill.attachmente       -0.922859   0.055651 -16.583  < 2e-16 ***
gill.attachmentf        0.816396   0.136750   5.970 2.37e-09 ***
gill.attachmentp       -2.450165   0.060717 -40.354  < 2e-16 ***
gill.attachments        0.142043   0.052091   2.727 0.006395 ** 
gill.attachmentx        0.089014   0.043392   2.051 0.040231 *  
gill.spacingd          -0.425196   0.036863 -11.534  < 2e-16 ***
gill.spacingf                 NA         NA      NA       NA    
gill.colore             2.304077   0.156835  14.691  < 2e-16 ***
gill.colorf                   NA         NA      NA       NA    
gill.colorg             0.833125   0.124899   6.670 2.55e-11 ***
gill.colork             0.854961   0.136405   6.268 3.66e-10 ***
gill.colorn             1.350254   0.119303  11.318  < 2e-16 ***
gill.coloro             1.131719   0.123598   9.156  < 2e-16 ***
gill.colorp             0.848246   0.122841   6.905 5.01e-12 ***
gill.colorr             0.896493   0.144108   6.221 4.94e-10 ***
gill.coloru             1.319431   0.145603   9.062  < 2e-16 ***
gill.colorw             0.819802   0.115149   7.120 1.08e-12 ***
gill.colory             1.818426   0.119373  15.233  < 2e-16 ***
stem.height             3.580359   0.196103  18.258  < 2e-16 ***
stem.width             -0.455339   0.225084  -2.023 0.043076 *  
stem.colore            17.897727 310.494233   0.058 0.954033    
stem.colorf            35.399323 334.067154   0.106 0.915610    
stem.colorg            15.906809 310.494228   0.051 0.959142    
stem.colork            19.545176 310.494262   0.063 0.949807    
stem.colorl            15.591779 310.494286   0.050 0.959950    
stem.colorn            17.491009 310.494224   0.056 0.955077    
stem.coloro            16.096825 310.494231   0.052 0.958654    
stem.colorp            19.131727 310.494251   0.062 0.950868    
stem.colorr            17.537012 310.494253   0.056 0.954959    
stem.coloru            17.419638 310.494223   0.056 0.955260    
stem.colorw            16.082601 310.494224   0.052 0.958691    
stem.colory            17.437089 310.494226   0.056 0.955215    
has.ringt              -0.006793   0.050561  -0.134 0.893124    
ring.typef             -0.911592   0.084553 -10.781  < 2e-16 ***
ring.typeg             -0.499679   0.101534  -4.921 8.60e-07 ***
ring.typel             -0.167323   0.100958  -1.657 0.097449 .  
ring.typem            -20.187930 230.159347  -0.088 0.930105    
ring.typep              0.894990   0.103503   8.647  < 2e-16 ***
ring.typer             -0.492646   0.105898  -4.652 3.29e-06 ***
ring.typez             16.741684  78.753199   0.213 0.831651    
habitatg                0.548824   0.041040  13.373  < 2e-16 ***
habitath                0.199731   0.065306   3.058 0.002225 ** 
habitatl               -0.501569   0.054900  -9.136  < 2e-16 ***
habitatm                0.103008   0.067655   1.523 0.127870    
habitatp               17.002090 230.326776   0.074 0.941156    
habitatu              -17.054326 378.536784  -0.045 0.964065    
habitatw              -17.436003 215.700665  -0.081 0.935574    
seasons                -1.618439   0.069225 -23.379  < 2e-16 ***
seasonu                 0.181688   0.025335   7.171 7.42e-13 ***
seasonw                -1.277528   0.049110 -26.014  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 67137  on 48854  degrees of freedom
Residual deviance: 45160  on 48776  degrees of freedom
AIC: 45318

Number of Fisher Scoring iterations: 16

La función aplicada, glm(), obtiene los valores residuales del modelo y los coeficientes de ajuste para cada una de las variables independientes. Además, se obtiene el p-value correspondiente para cada una de ellas. En el resumen mostrado mediante la función summary(), se pueden observar variables con dos o tres asteriscos, las cuales aportan bastante relevancia al modelo como predictores; sin embargo, las variables que poseen un solo asterisco o incluso ninguno, significa que apenas aportan relevancia a los resultados. A continuación, procedemos a predecir las clases del conjunto de entrenamiento y de validación. Tomando como umbral ‘0,5’, dividiremos los champiñones en comestibles y venenosos, de forma que si la probabilidad queda por encima de dicho umbral significará que el champiñón es comestible, y de lo contrario, si el resultado se mantiene por debajo, significará que el champiñón es venenoso. Para acabar, formaremos la matriz de confusión con los resultados de las predicciones para proceder a su análisis y valoración.

pred_train <- predict(rl_classiffier, newdata = training_set, type = "response")

Warning: prediction from a rank-deficient fit may be misleading

pred_train <- ifelse(pred_train > 0.5, "p", "e")
pred_train <- factor(pred_train, levels = c("e", "p"), labels = c("e", "p"))

confusion_m <- table(training_set$class, pred_train)
print(confusion_m)

   pred_train
        e     p
  e 16489  5256
  p  5573 21537

accuracy <- sum(diag(confusion_m)) / sum(confusion_m)
print(accuracy)

[1] 0.7783441

Una vez obtenido el resultado mostrado mediante la matriz de confusión, se observa que la predicción posee una precisión del 77,8%. A continuación, se muestra el resultado de aplicar el mismo proceso anterior sobre el conjunto de prueba.

pred_test <- predict(rl_classiffier, newdata = test_set, type = "response")

Warning: prediction from a rank-deficient fit may be misleading

pred_test <- ifelse(pred_test > 0.5, "p", "e")
pred_test <- factor(pred_test, levels = c("e", "p"), labels = c("e", "p"))

confusion_m <- table(test_set$class, pred_test)
print(confusion_m)

   pred_test
       e    p
  e 4087 1349
  p 1420 5358

accuracy_rl <- sum(diag(confusion_m)) / sum(confusion_m)
print(accuracy_rl)

[1] 0.7732929

Una vez obtenido este segundo resultado, podemos observar que hay una gran similitud entre los resultados obtenidos para ambos conjuntos de datos, resultando esta vez en un 77,3%. Antes de acabar de aplicar el clasificador de regresión logística, vamos a proceder a graficar la curva ROC, la cual nos aporta mayor visualización de la relación entre los falsos y verdaderos positivos.

library(ROCR)
pred_rl_roc <- prediction(as.numeric(pred_test), as.numeric(test_set$class))
perf_rl_roc <- performance(pred_rl_roc, "tpr", "fpr")
perf_rl_auc <- performance(pred_rl_roc, "auc")

print(perf_rl_auc@y.values[[1]])

[1] 0.7711691

plot(perf_rl_roc, col = "lightblue", lwd = 5)  

Observando la curva ROC resultante podemos comentar que se mantiene por encima de la diagonal, lo que es buena señal, pero se aproxima a ella, pudiendo haber proporcionado resultados mejores resulta ser un modelo bastante generalizado.

k-NN

Mediante el algoritmo de clasificación llamado k-NN se asigna una nueva observación a la clase más común entre sus “k” vecinos más cercanos en el espacio de características. Antes de nada, para poder aplicar el algoritmo k-NN, debemos transformar las variables categóricas en numéricas.

mushroom_num <- dummyVars(" ~ .", data = mushroom, fullRank = TRUE) %>% predict(mushroom)
mushroom_num <- as.data.frame(mushroom_num)

A continuación, se visualizan las variables categóricas ya codificadas y las dimensiones del dataset.

dim_mushroom <- dim(mushroom_num)
print(dim_mushroom)

[1] 61069    81

str(mushroom_num)

'data.frame':   61069 obs. of  81 variables:
 $ class.p              : num  1 1 1 1 1 1 1 1 1 1 ...
 $ cap.diameter         : num  0.24 0.262 0.221 0.223 0.23 ...
 $ cap.shapec           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.shapef           : num  0 0 0 1 0 0 1 0 1 1 ...
 $ cap.shapeo           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.shapep           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.shapes           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.shapex           : num  1 1 1 0 1 1 0 1 0 0 ...
 $ cap.surfacee         : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.surfaceg         : num  1 1 1 0 0 1 0 0 1 1 ...
 $ cap.surfaceh         : num  0 0 0 1 1 0 1 1 0 0 ...
 $ cap.surfacei         : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.surfacek         : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.surfacel         : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.surfaces         : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.surfacet         : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.surfacew         : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.surfacey         : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.colore           : num  0 0 0 1 0 0 0 1 0 1 ...
 $ cap.colorg           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.colork           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.colorl           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.colorn           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.coloro           : num  1 1 1 0 1 1 1 0 1 0 ...
 $ cap.colorp           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.colorr           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.coloru           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.colorw           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.colory           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ does.bruise.or.bleedt: num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.attachmentd     : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.attachmente     : num  1 1 1 1 1 1 1 1 1 1 ...
 $ gill.attachmentf     : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.attachmentp     : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.attachments     : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.attachmentx     : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.spacingd        : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.spacingf        : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.colore          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.colorf          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.colorg          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.colork          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.colorn          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.coloro          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.colorp          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.colorr          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.coloru          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.colorw          : num  1 1 1 1 1 1 1 1 1 1 ...
 $ gill.colory          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ stem.height          : num  0.5 0.53 0.525 0.465 0.487 ...
 $ stem.width           : num  0.164 0.175 0.171 0.154 0.166 ...
 $ stem.colore          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ stem.colorf          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ stem.colorg          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ stem.colork          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ stem.colorl          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ stem.colorn          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ stem.coloro          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ stem.colorp          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ stem.colorr          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ stem.coloru          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ stem.colorw          : num  1 1 1 1 1 1 1 1 1 1 ...
 $ stem.colory          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ has.ringt            : num  1 1 1 1 1 1 1 1 1 1 ...
 $ ring.typef           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ ring.typeg           : num  1 1 1 0 0 0 1 0 0 0 ...
 $ ring.typel           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ ring.typem           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ ring.typep           : num  0 0 0 1 1 1 0 1 1 1 ...
 $ ring.typer           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ ring.typez           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ habitatg             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ habitath             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ habitatl             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ habitatm             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ habitatp             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ habitatu             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ habitatw             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ seasons              : num  0 0 0 0 0 0 0 0 0 0 ...
 $ seasonu              : num  0 1 0 0 0 1 0 1 0 0 ...
 $ seasonw              : num  1 0 1 1 1 0 1 0 0 1 ...

Debido al hecho de necesitar la transformación de las variables categóricas en numéricas, procedemos en este paso a dividir los datos en los conjuntos de entrenamiento y de validación. En este caso, el valor 0 corresponde con los champiñones comestibles y el valor 1 con los champiñones venenosos.

library(caTools)
set.seed(18)

split <- sample.split(mushroom_num$class, SplitRatio = 0.8)
training_set_num <- subset(mushroom_num, split == TRUE)
test_set_num <- subset(mushroom_num, split == FALSE)

table(training_set_num$class)

    0     1 
21745 27110 

table(test_set_num$class)

   0    1 
5436 6778 

El siguiente paso será determinar el valor óptimo de k antes de proceder a aplicar la función. Para ello, utilizaremos el resultado que nos proporciona el cálculo de la raíz cuadrada del número de observaciones del conjunto de entrenamiento.

nrows_class <- NROW(training_set_num) 
k <- sqrt(nrows_class)
k <- round(k)
k

[1] 221

Una vez hemos obtenido el valor de k, procedemos a realizar las predicciones. Para llevar a cabo la aplicación del clasificador llamado ‘el vecino más cercano’, haremos uso de la función knn() de la librería class.

library(class)
set.seed(18)
pred_knn <- knn(train = training_set_num[, -1], test = test_set_num[, -1], cl = training_set_num$class, k = k)
summary(pred_knn)

   0    1 
5288 6926 

Una vez obtenidos los resultados de las predicciones en este caso, podemos de igual forma construir la matriz de confusión.

confusion_m <- table(test_set_num$class, pred_knn)
confusion_m

   pred_knn
       0    1
  0 5233  203
  1   55 6723

accuracy_knn <- sum(diag(confusion_m)) / sum(confusion_m)
accuracy_knn

[1] 0.9788767

Como resultado final de la matriz de confusión, obtenemos una precisión resultante del 97%, la cual mejora respecto al algoritmo de clasificación anterior, aunque tendiendo a poseer mayor sobreajuste. De forma más visual, obtenemos a continuación la gráfica de la curva ROC y el cálculo del área bajo la misma.

library(ROCR)
pred_knn_roc <- prediction(as.numeric(pred_knn), as.numeric(test_set_num$class))
perf_knn_roc <- performance(pred_knn_roc, "tpr", "fpr")
perf_knn_auc <- performance(pred_knn_roc, "auc")

print(perf_knn_auc@y.values[[1]])

[1] 0.9772709

plot(perf_knn_roc, col = "lightblue", lwd = 5)  

Una vez obtenemos la curva ROC y su área, observamos que el resultado es muy positivo en cuanto a precisión de los resultados, ya que como se puede observar gráficamente se aleja de la diagonal.

Clasificación con Árbol de Decisión

Los árboles de decisión se basan en la construcción de reglas lógicas (divisiones de los datos entre rangos o condiciones) a partir de los datos de entrada. Para trabajar con este clasificador comenzamos aplicando la función del árbol de decisión, rpart(), sobre el conjunto de datos de entrenamiento como se muestra a continuación.

library(rpart)
set.seed(18)
dt_classiffier <- rpart(class ~ ., data = training_set)

Una vez obtenido el resultado, lo graficaremos para facilitar así el análisis del resultado.

library(rpart.plot)
rpart.plot(dt_classiffier)

Construimos la matriz de confusión con los resultados obtenidos anteriormente para este caso.

pred_dt <- predict(dt_classiffier, newdata = test_set, type = "class")

confusion_m <- table(test_set$class, pred_dt)
confusion_m

   pred_dt
       e    p
  e 4257 1179
  p  804 5974

accuracy_dt <- sum(diag(confusion_m)) / sum(confusion_m)
accuracy_dt

[1] 0.8376453

Una vez tenemos el resultado para el valor de precisión en la predicción aplicando el árbol de decisión, 83%, procedemos a construir la gráfica de la curva ROC y a calcular el valor correspondiente al AUC. En este caso, continúa siendo mejor resultado que el obtenido mediante la regresión logística, puesto que resulta en una mayor precisión, pero al igual que el algoritmo de k-NN, tiende a ser más sobreajustado.

library(ROCR)
pred_dt_roc <- prediction(as.numeric(pred_dt), as.numeric(test_set$class))
perf_dt_roc <- performance(pred_dt_roc, "tpr", "fpr")
perf_dt_auc <- performance(pred_dt_roc, "auc")

print(perf_dt_auc@y.values[[1]])

[1] 0.8322468

plot(perf_dt_roc, col = "lightblue", lwd = 5)  

Clasificador Random Forest

El algoritmo de Random Forest trabaja mediante la combinación de árboles predictores tal que cada árbol depende de los valores de un vector aleatorio. Comenzamos aplicando dicho clasificador mediante la función llamada de la misma forma; es decir, randomForest(). En esta función, el valor correspondiente al parámetro llamado ‘ntree’ indica la cantidad de árboles de decisión que formarán parte del clasificador.

library(randomForest)
set.seed(18)
rf_classiffier <- randomForest(class ~ ., data = training_set, ntree = 250)

Mediante su gráfica, vamos a proceder a comparar los errores en función del aumento del número de árboles; es decir, cuanto más vaya aumentando el número de árboles hasta un umbral determinado, menor cantidad de errores poseerá la predicción.

plot(rf_classiffier)

Calculamos las predicciones sobre el conjunto de datos de prueba y construimos la matriz de confusión.

pred_rf <- predict(rf_classiffier, newdata = test_set, type = "class")

confusion_m <- table(test_set$class, pred_rf)
confusion_m

   pred_rf
       e    p
  e 5431    5
  p    0 6778

accuracy_rf <- sum(diag(confusion_m)) / sum(confusion_m)
accuracy_rf

[1] 0.9995906

Una vez obtenido el valor de la precisión para este caso, definimos la curva ROC y procedemos a realizar el cálculo del área bajo la curva. Mediante este último paso se puede observar su tan alta precisión, la cual indica demasiado sobreajuste sin ser conveniente.

library(ROCR)
pred_rf_roc <- prediction(as.numeric(pred_rf), as.numeric(test_set$class))
perf_rf_roc <- performance(pred_rf_roc, "tpr", "fpr")
perf_rf_auc <- performance(pred_rf_roc, "auc")

print(perf_rf_auc@y.values[[1]])

[1] 0.9995401

plot(perf_rf_roc, col = "lightblue", lwd = 5)  

Kernel SVM Classifier

El clasificador de la máquina vectorial, encuentra la curva que es capaz de separar y clasificar los datos de entrenamiento garantizando que la separación entre ésta y ciertas observaciones del conjunto de entrenamiento resulte ser lo mayor posible. Para llevar a cabo la aplicación del clasificador llamado ‘Máquina de Soporte Vectorial’ hacemos uso de la función svm(). En dicha función, los valores de los parámetros ‘type’ y ‘kernel’ hacen referencia al tipo de clasificador lo que significa que el kernel será de tipo radial y gaussiano.

library(e1071)
set.seed(18)
svm_classiffier <- svm(class ~ .,
  data = training_set,
  type = "C-classification", kernel = "radial"
)

A continuación, se calcula la predicción y se construye la matríz de confusión, las cuales resultan ser las siguientes.

pred_svm <- predict(svm_classiffier, newdata = test_set, type = "class")

confusion_m <- table(test_set$class, pred_svm)
confusion_m

   pred_svm
       e    p
  e 5129  307
  p  240 6538

accuracy_svm <- sum(diag(confusion_m)) / sum(confusion_m)
accuracy_svm

[1] 0.9552153

Como se puede observar, el valor de la precisión en la predicción en este caso resulta ser del 95%, un resultado bueno que no muestra señales de sobreajuste. Para finalizar, construimos la curva ROC correspondiente en este caso y calculamos su área.

library(ROCR)
pred_svm_roc <- prediction(as.numeric(pred_svm), as.numeric(test_set$class))
perf_svm_roc <- performance(pred_svm_roc, "tpr", "fpr")
perf_svm_auc <- performance(pred_svm_roc, "auc")

print(perf_svm_auc@y.values[[1]])

[1] 0.954058

plot(perf_svm_roc, col = "lightblue", lwd = 5)  

Conclusiones

Una vez aplicados los cinco distintos métodos de clasificación sobre nuestro dataset llamado ‘mushroom’ tras haber realizado antes su preprocesamiento, podemos concluir diciendo que el clasificador de Random Forest es el que ha resultado poseer un mayor valor en la precisión de la predicción en la clasificación y, por lo tanto, un menor valor para el error de predicción. Sin embargo, al aplicar este algoritmo hemos obtenido un mayor sobreajuste, el cual no beneficia al modelo ya que se busca que los resultados obtenidos sean precisos, pero también generalizados para los datos. Por esta razón, consideramos que resulta más beneficioso sacrificar parte del valor de precisión, teniendo en cuenta algunos valores de falsos positivos y falsos negativos, como ocurre en el caso de los clasificadores de la Máquina de Soporte Vectorial o k-NN, aprovechando así su capacidad de mayor generalización. Sobre los gráficos que se muestran a continuación se pueden comparar los distintos niveles de precisión y AUC para cada uno de los diferentes clasificadores que han sido aplicados.

accuracy_comp <- matrix(c(accuracy_rl, accuracy_knn, accuracy_dt, accuracy_rf, accuracy_svm), ncol = 5)

barplot(accuracy_comp,
  main = "Accuracy Comparison",
  xlab = "Accuracy (%)",
  ylab = "Method",
  names.arg = c("RL", "K-NN", "DT", "RF", "SVM"),
  col = "#7fd6d9"
)

perf_auc <- matrix(c(perf_rl_auc@y.values[[1]], perf_knn_auc@y.values[[1]], perf_dt_auc@y.values[[1]], perf_rf_auc@y.values[[1]], perf_svm_auc@y.values[[1]]), ncol = 5)

barplot(perf_auc,
  main = "AUC Comparison",
  xlab = "AUC (%)",
  ylab = "Method",
  names.arg = c("RL", "K-NN", "DT", "RF", "SVM"),
  col = "#7fd6d9"
)

Análisis aprendizaje no supervisado

En este apartado se analizará el dataset a través de algoritmos de aprendizaje no supervisado. En concreto, se probarán los algoritmos k-means y clustering jerárquico. Para ambos algoritmos, se seguirá el siguiente esquema:

Para ambos algoritmos se llevará a cabo el siguiente proceso:

1. Se representará gráficamente la distribución inicial de los datos en el dataset.
2. Se determinará el número óptimo de clústeres a utilizar para dividir los datos de forma adecuada.
3. Se representará gráficamente la distribución de los datos en función del número de clústeres elegido.
4. Se calculará el promedio de cada una de las variables en el dataset para cada uno de los clústeres resultantes.
5. Finalmente, y gracias a tener información a priori de la clase de cada champiñón, se calculará el “accuracy” del algoritmo, midiendo cómo de bien se está realizando la tarea de dividir los datos en clústeres de forma adecuada.

Para poder trabajar con algoritmos no supervisados será necesario que las variables sean numéricas. Para ello, se eliminarán las variables categóricas del dataset.

numerical_columns <- mushroom[, numerical_features]

Antes de comenzar con los algoritmos no supervisados, representaremos de forma gráfica la distribución inicial de los datos a través de un diagrama de dispersión 3D, dónde cada punto representa un champiñón y las variables que se representan son el diámetro del sombrero, la altura del tallo y el ancho del tallo. Cada variable está normalizada entre 0 y 1.

df <- as.data.frame(numerical_columns)

plot_ly(df,
  x = ~cap.diameter, y = ~stem.height,
  z = ~stem.width
) %>%
  add_markers(size = 1.5)

K-means

El algoritmo k-means es un método de clustering que permite dividir un conjunto de datos en k grupos o clústeres de manera que los puntos dentro de un mismo clústeres sean similares entre sí y diferentes a los puntos de los demás clústeres. La función kmeans() de la librería clústeres una implementación del algoritmo k-means en R.

Para utilizar la función kmeans(), es necesario especificar el número de clústeres que se desean obtener, que se indica a través del parámetro “centers”.

Por otro lado, el parámetro “nstar” indica el número de veces que se desea realizar el proceso de clustering. Cada vez que se ejecuta el proceso, se utiliza un conjunto diferente de semillas iniciales para los centroides de los clústeres y se obtiene un resultado diferente. Al especificar un valor para “nstart” mayor que 1, se obtienen varios resultados diferentes y se selecciona el que minimiza la suma de cuadrados total. Por tanto, establecemos “nstart” a 20 para que se realice el proceso 20 veces y se obtenga un resultado más robusto.

Una vez que se ha ejecutado la función con un determinado valor de “centers”, se puede calcular la suma de cuadrados internos (within groups sum of squares) para ese valor de “centers”. La suma de cuadrados internos es una medida de la variabilidad de los datos dentro de cada cluster. Cuanto mayor sea la suma de cuadrados internos, más dispersos estarán los datos dentro del clústery, por tanto, menos homogéneo será el cluster.

Para determinar el número óptimo de clústeres, se puede utilizar el método del codo, que consiste en representar la suma de cuadrados internos en función del número de clústeres y seleccionar el número de clústeres en el que se produce un “codo” en la gráfica. Este “codo” suele corresponder al punto en el que la disminución de la suma de cuadrados internos se vuelve más lenta y, por tanto, a partir del cual no se obtienen mejoras significativas en la calidad del clustering.

wss_per_k <- 0
for (i in 1:10) {
  kmeans_aux <- kmeans(numerical_columns, center = i, nstar = 20)
  wss_per_k[i] <- kmeans_aux$tot.withinss
}

Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)

par(mfrow = c(1, 1))
plot(1:10, wss_per_k,
  type = "b",
  xlab = "Number of clusters",
  ylab = "WSS",
)

Como se puede observar en la gráfica anterior, la suma de cuadrados internos disminuye a medida que aumenta el número de clústeres. Sin embargo, a partir de 2 clústeres, la disminución de la suma de cuadrados internos es más pequeña. Por lo tanto, se decide hacer uso 2 clústeres. En este caso específico, tiene sentido utilizar 2 clústeres, ya que conocemos que el dataset es binario.

Una vez determinado el número de clústeres, generamos el modelo de k-means.

km_model <- kmeans(df, center = 2, nstar = 20)

Para poder visualizar los resultados, se añade una nueva columna al dataset con el número de clúster al que pertenece cada observación. Una vez añadida, se puede representar la distribución de los datos en función de los clústeres obtenidos.

df$cluster<- factor(km_model$cluster)

plot_ly(df,
  x = ~cap.diameter, y = ~stem.height,
  z = ~stem.width, color = ~cluster
) %>%
  add_markers(size = 1.5)

Warning: minimal value for n is 3, returning requested palette with 3 different levels
Warning: minimal value for n is 3, returning requested palette with 3 different levels
Warning: minimal value for n is 3, returning requested palette with 3 different levels
Warning: minimal value for n is 3, returning requested palette with 3 different levels

Se puede observar que los champiñones de menor tamaño (en diámetro, altura y anchura) pertenecen al clúster 2 y los de mayor tamaño pertenecen al clúster 1.

A continuación, calcularemos el valor promedio de las variables para cada clúster generado con el modelo de k-means. Para ello, utilizaremos la función group_by() de la librería dplyr para agrupar los datos por clúster y la función summarise() para calcular el valor promedio de cada variable. Es importante destacar que, en este caso, los datos están normalizados, por lo que el valor promedio de cada variable no tiene un significado real. Sin embargo, nos permite comparar los valores de cada variable para cada cluster.

grouped_mushroom <- df %>%
  group_by(cluster) %>%
  summarise(
    mean_cap_diameter = mean(cap.diameter),
    mean_stem_height = mean(stem.height),
    mean_stem_width = mean(stem.width)
  )

grouped_mushroom

Observamos que los valores del clúster 1 son mayores que los del clúster 2, lo que indica que los champiñones del clúster 1 son de mayor tamaño que los del clúster 2.

A partir de este momento, hemos decidido modificar el dataset actual debido a que para aplicar técnicas como “silhouette” o “dendrogram” (para el caso de clustering jerárquico) es necesario que el dataset sea de menor tamaño. La función createDataPartition() de la librería caret en R permite dividir un conjunto de datos en dos grupos, uno de entrenamiento y otro de validación, de manera estratificada, es decir, manteniendo la proporción de elementos de cada clase en ambos grupos. En este caso, se está especificando que se desea mantener únicamente el 1% de los datos iniciales para el análisis, lo que implica que se está utilizando createDataPartition() para reducir el tamaño del conjunto de datos en lugar de para dividirlo en dos grupos. Al utilizar createDataPartition() de esta manera, se mantiene la proporción de elementos de cada clase en el conjunto de datos reducido.

Antes de reducir el dataset, es necesario convertir las variables categóricas en numéricas a través de variables dummy. Para ello, utilizamos la función dummyVars() de la librería caret para crear un objeto de tipo dummyVars y la función predict() para crear un nuevo dataset con las variables dummy.

mushroom <- dummyVars(" ~ .", data = mushroom, fullRank = TRUE) %>% predict(mushroom)
mushroom <- as.data.frame(mushroom)

set.seed(42)
split <- createDataPartition(mushroom$class, p = 0.01)
smaller_df <- mushroom[split$Resample1, ]

Comprobamos que la proporción de datos de cada clase se mantiene al hacer la partición.

initial_class_prop <- table(mushroom$class) / nrow(mushroom)
smaller_class_prop <- table(smaller_df$class) / nrow(smaller_df)

print(initial_class_prop)

        0         1 
0.4450867 0.5549133 

print(smaller_class_prop)

        0         1 
0.4386252 0.5613748 

Mostramos de forma gráfica las nuevas proporciones de la variable dependiente “class” en el dataset reducido.

print(ggplot(smaller_df, aes_string(x = smaller_df$class)) +
  geom_bar(fill = "#7fd6d9") +
  geom_text(stat = "count", aes(label = scales::percent(..count.. / nrow(smaller_df)), vjust = -0.25)) +
  labs(x = i, y = "Percentage") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)))

Una vez reducido el dataset, vamos a eliminar las variables categóricas para poder aplicar técnicas de clustering.

smaller_df <- smaller_df[, numerical_features]
dim(smaller_df)

[1] 611   3

Tras ejecutar la celda anterior se puede observar que el dataset ha pasado a tener 611 observaciones y 3 variables.

Como nos encontramos ante un dataset “nuevo”, en primer lugar, visualizaremos la distribución inicial de los datos a través de una gráfica 3D.

plot_ly(smaller_df,
  x = ~cap.diameter, y = ~stem.height,
  z = ~stem.width
) %>%
  add_markers(size = 1.5)

A continuación, vamos a estudiar cuál sería el número óptimo de clústeres para el dataset reducido haciendo uso de la medida de bondad interna “silhouette”. Para ello, utilizaremos la función fviz_nbclust de factoextra.

La medida silhouette toma valores entre -1 y 1. Un valor cercano a 1 indica que el punto está bien asignado al clúster y los puntos del clúster son muy similares entre sí. Un valor cercano a 0 indica que el punto está en un “área gris” y no está claramente asignado a ninguno de los dos clústeres. Un valor cercano a -1 indica que el punto está mal asignado al clúster y sería más apropiado para otro clúster.

La función fviz_nbclust() de la librería factoextra en R permite visualizar la medida silhouette para diferentes valores de “k” (número de clústeres) y ayudar a determinar el número óptimo de clústeres. Al utilizar esta función, se puede obtener un gráfico en el que se representa la medida silhouette en función del número de clústeres y seleccionar el valor de “k” en el que se obtiene el mayor valor de silhouette.

fviz_nbclust(smaller_df, FUNcluster = kmeans, method = "silhouette")

Según la gráfica, podemos afirmar que el número óptimo de clústeres es 2, ya que es el valor de “k” que maximiza la medida silhouette. También podemos observar que el valor de silhouette para “k” = 3 es cercano al obtenido para “k” = 2. Esto puede ser debido a que nuestro dataset, pese a ser binario, cuenta con datos muy dispersos.

Por último, visualizamos la distribución de los datos en función de los clústeres obtenidos.

km_sm_model <- kmeans(smaller_df, center = 2, nstart = 20)
cluster <- factor(km_sm_model$cluster)

plot_ly(smaller_df,
  x = ~cap.diameter, y = ~stem.height,
  z = ~stem.width, color = ~cluster
) %>%
  add_markers(size = 1.5)

Warning: minimal value for n is 3, returning requested palette with 3 different levels
Warning: minimal value for n is 3, returning requested palette with 3 different levels
Warning: minimal value for n is 3, returning requested palette with 3 different levels
Warning: minimal value for n is 3, returning requested palette with 3 different levels

Por último, calculamos el valor promedio de las variables para cada clúster generado con el modelo de k-means para el dataset reducido.

grouped_sm_mushroom <- smaller_df %>%
  mutate(cluster = cluster) %>%
  group_by(cluster) %>%
  summarise(
    mean_cap_diameter = mean(cap.diameter),
    mean_stem_height = mean(stem.height),
    mean_stem_width = mean(stem.width)
  )
grouped_sm_mushroom

Clustering jerárquico

Clustering jerárquico es un tipo de algoritmo de clustering que se utiliza para dividir a un conjunto de datos en grupos (clústeres) de forma que los datos en el mismo clústersean similares entre sí. La función hclust es una función en R que se utiliza para realizar clustering jerárquico.

Antes de aplicar el algoritmo hclust, es necesario calcular las distancias entre los puntos del conjunto de datos, para ello, utilizaremos la función dist de R. La función dist() calcula la distancia entre los puntos del conjunto de datos y devuelve una matriz de distancias. Por defecto, la función dist() utiliza la distancia euclídea para calcular las distancias entre los puntos del conjunto de datos.

Cabe destacar, que la función hclust hace uso por defecto del cáculo de distancia entre clústeres basado en el método de “Complete”. Este método de cálculo de distancia entre clústeres se basa en la distancia entre los puntos más lejanos de cada cluster.

distance <- dist(smaller_df)
hc_model <- hclust(distance)

Representamos el dendrograma para visualizar la distribución de los datos en función de los clústeres obtenidos.

dend_modelo <- as.dendrogram(hc_model)
plot(dend_modelo, ylab = "Similarity")

Hasta ahora, hemos obtenido la jerarquía de los datos, pero lo que realmente nos interesa es la clasificación de los datos en función de los clústeres. Cortaremos el dendrograma en un punto que nos interese para obtener los clústeres. En este caso, y a modo de prueba, hemos decidido cortar el dendrograma en 90 para obtener una visualización del dendograma cortado.

cut <- 0.9

dend_modelo %>%
  color_branches(h = cut) %>%
  color_labels(h = cut) %>%
  plot(ylab = "Similarity")

Para obtener el número óptimo de clúster, haremos uso de la medida interna de bondad silhouette. Para ello, utilizaremos la función fviz_nbclust de factoextra al igual que con k-means.

fviz_nbclust(smaller_df, FUNcluster = hcut, method = "silhouette")

Comprobamos que en este caso, el número óptimo de clústeres podría ser 2 o 3, ya que el valor de silhouette es muy similar para ambos casos. En este caso, hemos decidido utilizar 2 clústeres para poder comparar posteriormente los resultados con los obtenidos con el algoritmo de k-means.

Para generar el modelo de clustering jerárquico, utilizaremos la función cutree de R. Esta función nos permite generar el modelo de clustering jerárquico en función del número de clústeres que queramos obtener.

Calculamos la agrupación del modelo en función del número de clústeres que hemos decidido utilizar, y calculamos el valor promedio de las variables para cada clústergenerado.

jq_cluster <- cutree(hc_model, k = 2)

grouped_mushroom <- smaller_df %>%
  mutate(cluster = jq_cluster) %>%
  group_by(cluster) %>%
  summarise_all(mean)
grouped_mushroom

Visualizamos la agrupación de los datos en función de los clústeres obtenidos a partir del modelo de clustering jerárquico.

jq_clúster<- factor(jq_cluster)

plot_ly(smaller_df,
  x = ~cap.diameter, y = ~stem.height,
  z = ~stem.width,
  color = ~jq_cluster
) %>%
  add_markers(size = 1.5)

Con el objetivo de comparar los resultados obtenidos en los dos algoritmos, vamos a calcular el rendimiento de cada uno de ellos, haciendo uso del accuracy como medida de bondad externa.

En primer lugar, calculamos el accuracy del modelo de k-means. Supondremos que la clase 1 es la clase “e” y la clase 2 es la clase “p”. Para ello, obtenemos las clases reales y las clases predichas, y calculamos el accuracy.

Primero necesitamos volver a obtener el dataset reducido para poder tener las clases reales.

smaller_df <- mushroom[split$Resample1, ]

real_classes <- ifelse(smaller_df$class == "e", 1, 2)
predicted_classes <- km_sm_model$cluster
predicted_classes <- as.numeric(predicted_classes)

accuracy <- sum(real_classes == predicted_classes) / length(real_classes)
print(accuracy)

[1] 0.700491

Hacemos lo mismo con el modelo de clustering jerárquico, pero en este caso, supondremos que la clase 1 es la clase “p” y la clase 2 es la clase “e”.

real_classes <- ifelse(smaller_df$class == "e", 2, 1)
predicted_classes <- as.numeric(jq_cluster)

accuracy <- sum(real_classes == predicted_classes) / length(real_classes)
print(accuracy)

[1] 0.9885434

Tras comparar los resultados obtenidos en los dos algoritmos, podemos afirmar que el modelo de clustering jerárquico ha obtenido un accuracy mayor para este dataset, obteniendo un accuracy del 98% frente al 70% obtenido por el modelo de k-means.

Pese a obtener un accuracy mayor con el modelo de clustering jerárquico, hay que tener en cuenta que el objetivo de aprendizaje no supervisado es la agrupación de los datos en función de sus características, y no la predicción de una variable objetivo. Además, el accuracy obtenido con el modelo de clustering jerárquico es muy alto, lo que puede deberse a que el dataset utilizado es muy reducido y no presenta mucha variabilidad entre las clases.